#ifndef BLOSUM_H
#define BLOSUM_H

#include <boost/numeric/ublas/matrix.hpp>
#include <vector>
#include <set>
#include <map>
#include <string>
#include <fstream>

using boost::numeric::ublas::matrix;
using namespace std;

class Blosum {
    vector<string> sequences_;
    int sequenceLenght_;
    matrix<int> blosum_;
    set<char> alphabetSet_;
    map<char, int> alphabetMap_;
    int alphabetLen_;

public:
    /**
     * Konstruktor domyślny
     */
    Blosum ();
    
    /**
     * Dodawanie zbioru liter z którego składają się sekwencje
     * @param alphabet alfabet, na którym zbudowane są sekwencje uczące macierzy
     */
	  void setAlphabet(string alphabet );

    /**
     * Dodawanie zbioru sekwencji uczących macierz
     * @param file_path ścieżka do pliku z sekwencjami
     */
    void fillSequencesFromFile ( const char* file_path );

    /**
     * Dodawanie zbioru sekwencji uczących macierz
     * @param sequences sekwencje
     */
    void setSequences ( const vector< string >& sequences );

    /**
     * Obliczanie macierzy
     */
    void calculate();

    /**
     * Zwraca macierz BLOSUM
     */
    matrix< int > getBlosumMatrix() const;

    /**
     * Zapisuje macierz BLOSUM do pliku
     * @param ścieżka do pliku
     */
    void saveBlosumMatrix ( string file_path ) const;

private:
    /**
     * Zlicza ilość wystąpień każdego ze znaku (A, G, C, T) na pozycjach w sekwencjach.
     * Innymi słowy oblicza wektor n elementowy (n - długość pojedynczej sekwencji), gdzie
     * i-ty element to 4 liczby mówiące o tym, ile razy wystąpił na i-tej pozycji znak
     * odpowiednio (A, G, C, T).
     * Na przykład:
     * gdy 3-ci element wektora zawiera liczby {3, 4, 0, 1}, to znaczy, że w naszych sekwencjach
     * (pole sequences_) na 3-ciej pozycji znak A wystąpił 3 razy, znak G 4 razy, znak C 0 razy,
     * a znak T 1 raz.
     *
     * @param counts miejsce na umieszczenie wyników obliczeń
     */
    void calculateCounts ( vector< vector<int> > & counts ) const;

    /**
     * Oblicza macierz wystąpień elementów jak na slajdzie 33/48 z wykładu,
     * ale tylko jednej połowy, bo macierz w oryginale jest symetryczna
     * @param counts ilości wystąpień elementów na poszczególnych pozycjach
     * @param matrix1 miejsce na wynik obliczeń
     * @return suma elementów macierzy (wszystkich, nie tylko jednej połowy)
     */
    int calculateStep1Matrix ( matrix< int >& matrix1,
                               const vector< vector< int > >& counts ) const;

    /**
     * Oblicza macierz prawdopodobieństw jak na slajdzie 34/48
     * @param probabilities ta macierz zostanie wypełniona obliczonymi prawdopodobieństwami
     * @param marginals ten wektor wypełniony będzie prawdopodobieństwami brzegowymi
     * @param step1_matrix wyliczona w poprzednim kroku macierz z występowaniami znaków w sekwencjach
     * @param sum liczba wystąpień wszystkich znaków w sekwencjach
     */
    void calculateProbabilitiesMatrix ( matrix< double > & probabilities,
                                        vector<double> & marginals,
                                        const matrix< int > & step1_matrix,
                                        int sum ) const;

    /**
     * Oblicza macierz BLOSUM
     * @param probabilities_matrix macierz prawdopodobieństw
     * @param marginal_probabilities wektor prawdopodobieństw brzegowych
     */
    void calculateBlosum ( const matrix< double >& probabilities_matrix,
                           const vector< double >& marginal_probabilities );
    
    // utils
    inline int charToIndex ( char c ) const {
        return alphabetMap_.at(c);
    }

    inline char indexToChar ( int index ) const {
        for(map<char, int>::const_iterator i = alphabetMap_.begin(); i != alphabetMap_.end(); ++i) {
            if(i->second == index) {
                return i->first;
            }
        }
        return '?';
    }

    inline bool isInAlphabet ( char c ) const {
        return alphabetSet_.find ( toupper ( c ) ) != alphabetSet_.end();
    }

    // do testów
    void testShowCounts ( const vector< vector<int> > & counts ) const;
    void testShowStep1 ( const matrix< int > & substitutions, int sum ) const;
    void testShowProbabilities ( const matrix< double >& m, vector<double> & marginals ) const;
    void testShowBlosum();
};

#endif // BLOSUM_H
